% 主要参数设置
g = 9.81;             % 重力加速度 (m/s^2)
rho = 1.225;          % 空气密度 (kg/m^3)
Cd = 0.8;             % 空气阻力系数 (假定)
A = [0.5, 0.4];  % 各运动员迎风面积 (m^2)
m = [60, 70 ];    % 各运动员质量 (kg)
t_start = [0, 0.3 ]; % 各运动员起跳时间 (s)
k = [5000, 5200];  % 站立位置对应弹簧刚度 (N/m)
alpha = 1e5;          % 非线性项系数 (N/m^3)
dt = 0.001;           % 计算步长 (s)
T_end = 20;           % 仿真总时间 (s)
time = 0:dt:T_end;    % 时间轴
% 预分配变量
num = length(m);
z = zeros(num, length(time));      % 各运动员位移
v = zeros(num, length(time));      % 各运动员速度
F_spring = zeros(1, length(time));% 蹦床总弹力
% 数值积分求解动力学
for i = 1:num
    % 每个运动员分别计算
    idx0 = find(time >= t_start(i), 1);
    % 初始条件
    z(i, idx0) = 0;  v(i, idx0) = 0; 
    for t = idx0:length(time)-1
        % 当前速度和位移
        vi = v(i,t); zi = z(i,t);
        % 空气阻力 (取与速度反向)
        Fair = 0.5 * rho * Cd * A(i) * vi^2 * sign(-vi);
        % 弹簧力
        if zi < 0
            Fs = -( k(i)*zi + alpha*zi^3 );
        else
            Fs = 0;
        end
        % 运动员动力学方程 (欧拉积分示例)
        a = -g - Fair/m(i) + Fs/m(i);
        v(i,t+1) = v(i,t) + a*dt;
        z(i,t+1) = z(i,t) + v(i,t)*dt;
    end
end
% 叠加总弹簧力
for t = 1:length(time)
    F_spring(t) = sum( k .* min(z(:,t),0) ); % yi<0才贡献力
end
